Network-constrained 1st-Order Spatial Point Pattern Analysis

Author

Tang Xin Yi

Packages

Code
pacman::p_load(sf,tidyverse,spNetwork,tmap,classInt, viridis)

Import data

Prototype was done focusing on provinces of Thailand. Yet in the actual shiny app due to technical restrictions, we changed to be targeting at Bangkok districts only.

Code
crime_data <- read_rds("../data/rds/accidents_thai.rds") %>% 
  st_transform(crs = 32648)
thai_bound <- st_read(dsn = "../data/thai_adm_boundary",
                      layer = "tha_admbnda_adm2_rtsd_20220121") %>% 
  st_transform(crs = 32648)
Reading layer `tha_admbnda_adm2_rtsd_20220121' from data source 
  `/Users/tangtang/Desktop/IS415 Geospatial Analytics and Applications/project/gga-ThaiRoad/data/thai_adm_boundary' 
  using driver `ESRI Shapefile'
Simple feature collection with 928 features and 19 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 97.34336 ymin: 5.613038 xmax: 105.637 ymax: 20.46507
Geodetic CRS:  WGS 84
Code
thai_road <- read_rds("../data/rds/thai_roads.rds") %>% 
  st_transform(crs = 32648)

Visualising difference in network constrained kernel density through different filters by columns

Some potential filtering to apply to dataset to allow user to zoom in into specific area he/she is interested in:

  • Accident categories

  • Fatal accident [yes/no]

  • Province (MUST-HAVE)

But let’s visualised the overall look

Code
plot(crime_data)

Code
table(crime_data$fatal_accident)

  no  yes 
5974  107 
Code
table(crime_data$accident_categories)

    driver_factors   external_factors             others           speeding 
               413                167                506               4857 
traffic_violations 
               138 
Code
table(crime_data$province_en)
< table of extent 0 >

Example: filter by province = “Bangkok”, accident_categories = “speeding”, fatal_accident == “yes” (Allowing users to change freely)

Code
province_factor = "Bangkok"
accident_category = "speeding"
fatal = "yes"

crime_data_filtered = crime_data %>% 
  filter(accident_categories == accident_category,
         fatal_accident == fatal)

Network-constrained Kernel Density Estimation

Get the portion of road network for the selected province

Code
thai_prov <- thai_bound %>% 
  filter(ADM1_EN == province_factor)
Code
tmap_mode("view")
tm_shape(thai_prov) +
  tm_polygons() +
tm_shape(crime_data_filtered) +
  tm_dots() +
tm_shape(thai_road) +
  tm_lines()
Code
thai_road <- thai_road %>%
  filter(fclass %in% c("motorway","primary","secondary","tertiary"))
road_prov <- st_intersection(thai_road, thai_prov)
road_prov_sf <- st_cast(st_collection_extract(road_prov,
                                             type = "LINESTRING"),"LINESTRING")

Prepare lixel object

1st try with 1000m diameter

Code
lixels <- lixelize_lines(road_prov_sf,
                         1000,
                         mindist = 500)

Generate line centre points

Code
samples <- lines_center(lixels)
densities <- nkde(road_prov_sf,
                  events = crime_data_filtered,
                  w = rep(1, nrow(crime_data_filtered)),
                  samples = samples,
                  ## important
                  kernel_name = "quartic",
                  bw = 500,
                  ## important
                  div = "bw",
                  method = "simple",
                  digits = 1,
                  tol = 1,
                  grid_shape = c(1,1),
                  max_depth = 8,
                  # aggregate events within 10m radius (faster calculation)
                  agg = 10,
                  sparse = TRUE,
                  verbose = FALSE
                  )

Generate plot

Code
samples$density <- densities
lixels$density <- densities

samples$density <- samples$density*1000
lixels$density <- lixels$density*1000
Code
tmap_mode('view')
tm_shape(lixels) +
  tm_lines(col="density",
           lwd = 3)
Code
tmap_mode('plot')

UI design for Analysis

Storyboard of NKDE

The user would be able to filter the dataset to see impact different categories on the network-constrained kernel density map created and use the interpretation rules given to derive at a conclusion.